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Abstract 

We show that the plastic deformation of snow under uniaxial compression is characterized by complex 
spatio-temporal strain localization phenomena. Deformation is characterized by repeated nucleation 
and propagation of compaction bands. Compaction bands are also observed during the very first stage 
of compression of solid foams where a single band moves across the sample at approximately constant 
stress. However, snow differs from these materials as repeated nucleation and propagation of bands occurs 
throughout the subsequent hardening stage until the end of the deformation experiment. Band nucleation 
and/or reflection of bands at the sample boundaries are accompanied by stress drops which punctuate 
the stress strain curve. A constitutive model is proposed which quantitatively reproduces all features 
of this oscillatory deformation mode. To this end, a well-established compressive plasticity framework 
for solid foams is generalized to account for shear softening behavior, time dependence of microstructure 
(‘rapid sintering’) and non-locality of damage processes in snow. 


1 Introduction 

Irreversible deformation of snow plays an important role in a number of problems ranging from the interaction 
of snow with winter sports equipment [1] over vehicle traction to snowpack stability and avalanche release 
[2]. The strength of snow is strongly influenced by the presence of bonds connecting ice granules in the 
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snow microstructure. As a consequence, shear deformation reduces strength as bonds are broken, but 
densification increases strength as higher density leads to an increased number of contacts. The formation 
of bonds between adjacent ice granules can be envisaged as a thermodynamically driven sintering process [3] 
which leads to strengthening of snow over time (ageing) [4]. The interplay between strain softening and age 
hardening leads to ductile behavior at low strain rates, when broken bonds have sufficient time to re-form, 
and to quasi-brittle behavior at high strain rates when this is not possible. This transition can be observed 
both in pure shear [4] and in compressive deformation [5]. 

At intermediate strain rates where ageing and softening processes occur on the same time scale, complex 
spatio-temporal deformation patterns and oscillation phenomena may occur. Such phenomena are well 
documented in deformation of metals where they have been analysed under the generic term of strain rate 
softening instabilities [6, 7, 8] but have never been studied in snow or similar cohesive-granular materials. 
The present investigation presents the first experimental and computational investigation of spatio-temporal 
deformation oscillations in snow as observed in compression experiments which are carried out at strain rates 
that are borderline between the ductile and quasi-brittle deformation regimes. 

2 Experimental method 

Specimens of artificially produced dry snow with density p = 370kg/nW 3 and mean grain size £ ss 0.2mm 
were contained within a rectangular transparent container with Aluminium alloy side walls and a glass front 
and back. The specimens were compacted from above by an anvil moving at fixed rates of 1.125 mms' 1 and 
5mms _1 , thus providing nominal strain rates of e ex t = 5.77 x 10 _3 s _1 and 2.56 x 10 -2 s -1 . The experiments 
were carried out at a temperature of T — — 10°C. During deformation, the driving force was recorded by 
a load cell located above the anvil. A 18 megapixel camera was located in front of the specimens with 
illumination provided by a flashgun and diffuser located behind the specimens. Images of the transmitted 
light were recorded at 0.25 s intervals. When the images are viewed, compacted areas become apparent as 
darker regions (Fig. 1, left). It was generally found that compaction occurs in a heterogeneous manner by 
motion of compaction fronts which separate regions of different density (Fig. 1, see also results section). To 
quantify this phenomenon in terms of displacements, each image series was analyzed using digital image 
correlation (DIC) software. Local strain and strain rate tensors were then calculated from the displacement 
fields. Further details of the experimental set-up are given in the Supplementary Material. 
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Figure 1: Strain and strain rate pattern in a specimen deformed at e ex t = 5.77 x 10 _3 s _1 to 9.94 % overall 
compressive strain; left: photographic image showing the compaction front, the arrow indicates the direction 
of front propagation; center: corresponding strain distribution as obtained by DIC; right: corresponding 
strain rate distribution; note that the outermost edges of the sample are not covered by the DIC analysis. 

3 Constitutive Model 


To model compressive deformation of snow, we start out from a model proposed by Zaiser et al. [9] who 
generalize the compressible plasticity model formulated by Deshpande and Fleck [10] for deformation of 
solid foams to include non-locality and strain localization phenomena. We modify this model to account for 
the interplay between shear softening of snow associated with the breaking of bonds between ice granules 
upon shear deformation [11], and rapid sintering processes which re-establish intergranular links and restore 
strength over time. Pore pressure effects, which may influence dynamic snow compaction under impact loads 
[1], are not relevant at the presently investigated compaction rates and are therefore not taken into account. 

The stress tensor is denoted as cr with components a t j. It fulfils the momentum balance equation which 
for quasi-static deformation processes as considered here is = 0 where we sum over repeated indices. 
The first invariant of the stress tensor is Tr(cr) = an, and the deviatoric stress is (r^ev — cr — (l/3)Tr(<x)l 
with components <7^ ev . Our model is based on the Deshpande-Fleck yield criterion which has been shown to 
give a good description of the deformation of solid foams of low to intermediate density [10] and may thus 
be appropriate for snow where typical relative densities are in the range of 0.1-0.4. It is given by 
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where <r eq is the equivalent stress and the elastic domain is defined by / < 0. The plastic strain rate 
is assumed to derive from an associated flow rule, e pl = e(df /dcr) where the equivalent strain rate e is 
determined by imposing that, under temporally changing constraints, the system must not enter the domain 
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/ > 0 . The non-dimensional parameter a in Eq. 1 characterizes the degree of volumetric change during 
deformation. For the limit case a — 0, the model reduces to the standard Von Mises plasticity thory. 

We consider deformation at low to intermediate strains and use an additive decomposition of the strain 
tensor e = e pl + e el into elastic (reversible) and plastic (irreversible) contributions. The elastic strain e el 
is related to the stress via Hooke’s law where we assume isotropic elastic behavior with Young’s modulus 
E and Poisson number v. The elastic modulus E in general depends on density and may increase during 
densification. To describe this dependency, we assume a power-law density dependence as typically observed 
for solid foams [12], 

E = E 0 (p/p 0 ) c , (2) 

where p is the stress-free density of the snow, po is a reference density which we take to be the density prior 
to compaction, and Eq is the corresponding Young’s modulus. 

To complete our constitutive framework we need to relate the flow stress S to the plastic strain tensor. 
We assume S to depend both on the density of the snow and on its structure. Phenomenologically, also the 
yield stress of many solid foams has a power-law density dependence [12]: 

5 - E(p/po)°, (3) 

where the pre-factor S depends on the snow microstructure and may be affected by structural softening 
under shear deformations. The evolution of the stress-free density p as a function of plastic strain follows 
from the continuity equation dp/dt + pTr(e) pl = 0 which can be integrated to p — po exp[— £y] where the 
volumetric plastic strain is given by £y — Tre pl . Combining with Eq. 3 gives 

S- fc Eexp[— a£y], (4) 

which for negative volumetric strain describes densification-induced strengthening. 

Shear deformation of snow is associated with failure of bonds connecting ice granules and thus reduces 
the load bearing capacity of the snow microstructure. We note in passing the proposal of Heierli and 
Zaiser [13] and Heierli et. al. [14] that volumetric compaction may also be associated with (transient) 
structural softening. However, the simple viewpoint that densification implies strengthening while shear 
implies softening, which we adopt in the following, turns out to be fully adequate for understanding the 
experiments. We characterize the magnitude of shear deformation in terms of the equivalent shear strain 
£g* which relates to the second invariant of the deviatoric part of the plastic strain tensor by (eg 1 ) 2 — 
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(2/3)e?j’ dev e? 1,dev . The changes in strength which result from the competing processes of shear-induced 
bond failure and strength recovery due to rapid sintering processes are described by the phenomenological 
equations 

i pl 1 

£ = £ r + £ s , d t S s = -^S s + -(S°-£ s ). (5) 

£S T 

In these equations, Sr is the residual strength associated with intergranular friction which remains even 
when all bonds between ice granules have failed, £g is the structural strength contribution of intergranular 
bonds, £g is a characteristic shear strain associated with bond failure, £g is the limit strength reached 
after prolonged ageing, and r is the characteristic time constant for strength recovery by sintering/ageing 
processes. 

Owing to the possibility of structural softening, deformation may proceed in a localized manner even if the 
external loading is homogeneous, and spontaneous strain localization is indeed manifest in our experimental 
data. In such situations, the mathematical formulation of the deformation problem may become ill-posed, 
leading to mesh dependence of numerical solutions. To mitigate this problem, we adopt the suggestion of 
Aifantis [15] of adding a second-order gradient of the plastic shear strain to the yield function, a method 
which has been proven to be thermodynamically consistent by Gurtin and Anand [16] (for alternative 
regularization methods, see [17]). This leads to the constitutive equation 

S = (Sr + £ s ) exp[—a£y] + EofAef . (6) 

Mathematically, the non-local term suppresses strain localization on scales below the characteristic length 
l and defines a finite, mesh independent width of the deformation bands. The Laplacian in the yield 
condition necessitates a corresponding higher-order boundary condition, which we take to be ro.VEg 1 = 0, 
i.e., we require the gradient of the shear strain to be zero in the direction n perpendicular to the specimen 
surface. This choice of the higher-order boundary condition is found to correctly reproduce the behavior 
observed when a deformation band reaches the end of the specimen. 

In applying our model to the compression experiments considered here, we assume the parameters 
v — 0, a = 1/2. This allows us to reduce the model to three coupled equations for the axial compressive 
stress s a xx where x is the coordinate along the compression axis, the axial strain e := £ xx < 0, and for 
the structural strength Sg. A derivation of the simplified equations and a demonstration of the correctness 
of the assumptions ;/ ss 0, a ss 1/2 is given in the Supplementary Material. The axial stress is found to be 

s = E eS (e ext -(e)), with E~ j 1 = (E~ 1 (e(x))) (7) 
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where E e s is the effective elastic modulus of the sample. (...) denotes the spatial average along the x axis. 
The external strain e ex t, which is negative, is changed at a constant rate, and the concomitant change in 
the plastic strain e follows from the yield condition, 


s < 


(E r + Eg) exp[—ae] + E 0 l 2 


d 2 e 
dx 2 


( 8 ) 


where Eg evolves according to 

<9tE s = —Eg + - (Eg — Eg). (9) 

es t 

From these equations, a necessary condition for structural softening and strain localization derives as 
aeg < Eg/(ER + Eg). This condition (for derivation and further discussion, see supplementary material) is 
tantamount to the requirement that the strain softening processes due to shear deformation must outweigh 
the competing process of densification-induced hardening. 

To solve the system of integrodifferential equations for e, s and Eg, we start from initial conditions 
s — e = 0, Eg = Eg. In each time step we change e ex t according to the imposed strain rate and evaluate s 
and e self-consistently from Eqs. 7 and 8 with the boundary condition <9 x e(0) = d x e(L) = 0. By comparing 
with the previous time step we determine the local plastic strain rates e, evaluate the changes in Eg from 
Eq. 9, and repeat until a given end strain is reached. 


4 Results 

4.1 Experimental Observations 

Fig. 1 shows a snapshot of the patterns of axial strain e xx (x,y) and strain rate e xx (x,y) emerging during a 
compression test. The plots use Eulerian coordinates, i.e., each volume element is shown at its position in 
the laboratory frame. It is seen that compaction proceeds in a strongly heterogeneous manner, as a region 
of increased density at the top of the sample is separated from an uncompacted region at the bottom by a 
sharp moving front where the strain rate is concentrated. This front is perpendicular to the compression 
direction. Owing to the intrinsic structural heterogeneity of the snow microstructure, the deformation state 
behind the front is not completely homogeneous. To facilitate comparison with our model, which provides 
an effective one-dimensional description, in the following we average the strain and strain rate fields over 
the y direction. 

To visualize band motion, we use space-time plots where we plot the y-averaged strain rate etot (x,t) = 
(txx{x,y,t)) y in colorscale as a function of x and t (Fig. 2). In these plots we use Lagrangian x coordinates. 
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Figure 2: Top: space-time plots of the evolution of deformation activity as deduced from DIC; bottom: space- 
time plots of deformation activity as obtained from simulation; center: corresponding force vs. time/strain 
curves (black: experiment, red: simulation), the simulated curves have been corrected for friction effects as 
detailed in the Supplementary Material, (a) represents a sample deformed at e ex t = 5.77 x 10 -3 s -1 , (b,c) 
samples deformed at e ex t = 2.56 x 10~ 2 s _1 ; experimental data for deformation of samples (a) and (b) are 
also shown in Supplementary Movies 1 and 2. 
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i.e., x denotes the position along the compression axis in the initial configuration (thus the x interval occupied 
by the sample in the graphs does not shorten during compression). The colour contrast demonstrates the 
strongly localized nature of deformation. A moving deformation front appears on the space-time plot as an 
inclined zone of high strain rate. The inclination slope defines the front velocity in the Lagrangian frame. 
The first, front nucleates at the top of the sample and moves at constant (Lagrangian) speed downwards. 
Once it reaches the bottom of the sample, a new front nucleates there and moves upwards across the sample 
until it reaches the top. Repetition of this process leads to a bouncing motion of the locus of deformation. 
Oscillatory features are also manifest on the force vs. time curves as band nucleation is associated with an 
up-down oscillation of the driving force (see the dashed lines in Fig. 2 which illustrate the correspondence). 
Occasionally the bouncing pattern is interrupted as seen in Fig. 2(c) where the first band gets stuck in the 
middle of the sample and deformation is taken over by a second band nucleating at the bottom and moving 
upwards until it merges with the first, whereafter a new band nucleates at the top and the bouncing pattern 
is resumed. The observations are similar for both investigated deformation rates. 

4.2 Simulation results 

In our simulations we use the following model parameters (for discussion, see the Supplementary Material): 
Initial elastic modulus Eq — 10 MPa, volumetric change parameter a — 0.499, friction coefficient at container 
surface /i = 0.5, mean initial (homogeneous) strength of snow So = 8 x 10 4 Pa, exponent of density 
dependence of elastic modulus c == 3, exponent of density dependence of strength a — 8, softening strain 
£s = 0.05, residual strength Sr = 0.05Sq, characteristic time for rapid sintering r == 10s. The internal length 
is taken to be £ = 0.2mm, which we also use as spacing of our computational grid. To mimic microstructural 
heterogeneity, the initial strength at the different grid points is assumed to be Weibull distributed with 
mean Sq and Weibull exponent fi — 5. Finally, to mimick the effect of near-surface heterogeneities and 
stress concentrations due to uneven loading at the specimen ends, structural strength is assumed to decrease 
towards zero across two surface layers of width 1mm (5 grid points) and 0.4mm (2 grid points), located at 
the top and bottom of the specimen, respectively. 

Simulations were carried out at imposed strain rates e ex t — 5-7 x 10 _3 s _1 and e ex t = 2.5 x 10 _2 s _1 , 
matching the experimental data shown in Fig. 2. The two runs shown for the higher strain rate differ only 
with respect to the initial random distribution of local strength. 

It is seen that, similar to the experiment, the nucleation of new bands is associated with oscillations in 
the force vs. time/strain curves, though these are less pronounced than in the experimental counterparts. 
In the absence of strength fluctuations we consistently observe a bouncing band pattern similar to Fig. 2(a). 



Strength fluctuations occasionally induce band arrest (first dashed line in Fig. 2(b)) and/or intermittent 
band propagation, as also observed experimentally. Neither in the simulations nor in the experiments do we 
see a strong influence of the different band propagation modes on the overall shape of the force-displacement 
curves. 

Varying the strain rate by two orders of magnitude, we see a transition between qualitatively different 
types of spatio-temporal deformation patterns (Fig. 3) depending on the value of the product e ex t t. For 
e ex tT ^ 1, we find a single compaction band followed by homogeneous deformation, a behavior that is 
typical for solid foams [9]. For e ex i t — 1 we see first one secondary band and then, at even lower strain 
rates, a procession of bands similar to those observed in our experiments (see also the discussion of stability 
criteria in the Supplementary Material). 



Position (1 -x/h) 


Figure 3: Strain pattern evolution for different imposed strain rates, compressive strain profiles are plotted 
as a function of the position along the compression axis (the top of the sample is at the left of the graphs); 
adjacent profiles differ by equal total strain (time) increments, (a) e ex tr = 10, (b) e ex tT = 1, (c) e ex tT = 0.1, 
for other parameters see text. 


5 Discussion and Conclusions 

We have for the first time demonstrated that the irreversible deformation of snow may proceed in an 
intrinsically inhomogeneous manner characterized by the repeated nucleation and propagation of localized 
deformation bands. We could relate this phenomenon to the competition between structural softening 
under deformation, and age hardening due to rapid sintering processes which restore strength over time. 
Our findings clearly demonstrate the necessity to investigate and account for strain localization phenomena 
even when considering initially homogeneous snow samples under homogeneous loads. If shear softening 
outweights strengthening due to densification, unstable deformation accompanied by spatio-temporal strain 
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localization is bound to occur, otherwise, deformation is expected to proceed in a homogeneous manner 
associated with monotonic hardening. The difference is manifest not only in laboratory experiments, but in 
everyday experience: Walking on snow is, in the second case, associated with a foot resistance that increases 
with increasing penetration depth - as a consequence, we can walk on the snow. In the strain softening case, 
on the other hand, once the critical load needed to nucleate a compaction band is exceeded the foot travels 
downwards together with the compaction band at constant or even decreasing load - we have the feeling 
of breaking through the snow. This behavior, commonly interpreted in terms of snow heterogeneity ("there 
is a crust on the snow"), in fact reflects the intrinsic interplay of strain softening and strain localization 
characteristic even of homogeneous snow. 

The formation of compaction bands has close analogies in deformation of metallic foams, where the 
initial stage of compressive deformation is often characterised by localized deformation bands which can 
be explained in terms of localized structural softening (buckling or failure of cell walls) [17, 9]. However, 
repeated nucleation and propagation of deformation bands, which is a conspicuous feature of the present 
experiments, has never been observed in solid foams. The difference is readily understood if we consider 
the extremely high homologous temperatures in the present experiments where rapid sintering may lead to 
strength recovery on short time scales [18, 19]. As a consequence, if the characteristic time required for 
propagation of a deformation band is on the same order of magnitude as the characteristic time of strength 
recovery in a snow sample, repeated band nucleation may be observed. 

In the context of snow, the interplay between softening and ageing has been modelled by a number of 
authors (e.g. Louchet [20 \ 1 Reiweger et. al. [21]. However, models published to date cannot describe spatio- 
temporal strain localization, either because they implicitly assume homogeneous deformation [20] or because 
they consider spatial couplings within a mean-field framework which does not account for spatial structure 
of the deformation field [21]. We have demonstrated in the present paper that generalized continuum models 
(here: gradient plasticity), combined with a simple phenomenological description of strain softening and age 
hardening processes, offer a framework capable of quantitatively describing the observed spatio-temporal 
phenomena. Beyond the scope of the present study, such models may offer a generic framework to describe 
snow deformation also under heterogeneous and time-dependent loading conditions. 
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1 Experimental method 

Ice chips were cooled to — 80°C and crushed into a dry powdery material using a stan¬ 
dardised process [1], This material was then placed in a covered container in the cold lab 
for four hours to allow it to reach the temperature of — 10°C at which the deformation 
experiments were carried out. This resulted in a cohesive aggregate of ice granules with 
density p — 370 kg/m -3 and grain sizes ranging from 0.1 to about 0.5 mm (mean grain 
size £ ~ 0.2 mm). As is typical for laboratory- or machine-made snow, the density of our 
samples is at the upper end of the range observed for natural snow. The artificial snow 
thus produced was sifted through a 1.68 mm aperture sieve into a flat rectangular box 
with aluminium alloy side walls and a glass bottom. Excess snow was removed with a 
straight edge and the box was closed with a second sheet of glass screwed on top of the 
specimen. The container was then rotated upright to provide a specimen w — 250 mm 
wide, L — 195 mm high and d — 20 mm thick, enclosed in a container with transparent 
front and back walls within an apparatus as shown in Fig. 1. The total time from sieving 
to deformation was approximately 5 minutes. 

The specimen was compacted by an anvil moving downwards at fixed rate. During 
deformation, the driving force was recorded by a load cell located above the anvil. The 
recorded force data were filtered to remove high frequency oscillations associated with 
the motor drive. Typical compressive forces at the end of a deformation experiment 
(compressive strain e xx ~ 0.3) were of the order of 250 N. As the specimen is compacted, 
it exerts lateral forces on the side walls of the container. These forces were evaluated 
by determining, at the end of a deformation experiment, the sideways deflection of the 
container side walls using a pair of dial test indicators, and comparing with calculations 
made under the assumption that the response of the side walls is that of an elastic thin 
plate, with elastic properties determined by a calibration experiment. This yielded a 
total lateral force of about 45 N at compressive strain e xx = 0.3. 

An 18 megapixel camera was located in front of the specimen with illumination pro¬ 
vided by a flashgun and diffuser located behind the specimen. Images of the transmitted 
light were recorded at 0.25 s intervals. To analyze spatio-temporal deformation patterns, 
each image series was analysed using the digital image correlation software GeoPIV [2]. 
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This software takes samples of an image (known as patches) and searches for them in 
subsequent images, thus producing a displacement vector for each patch. The software 
is capable of locating a patch with sub-pixel accuracy. When this software was used to 
track small translations of an undeformed specimen, the mean error was found to be sig¬ 
nificantly less than 40 /am. The images were analyzed sequentially, using a grid of 64x64 
pixel patches spaced at 4.5 mm and with the patch positions reset to this grid for each 
image. Local strain tensors were then calculated from the displacement fields by taking 
the tensorial displacement gradient which was evaluated using the ParaView package [3], 
and strain rate tensors were evaluated from strain tensors pertaining to sequential images 
in terms of finite differences. 


2 Application of the constitutive model to the com¬ 
pression experiments 


The general constitutive framework has been defined in the main paper. When applying 
this framework to the compression experiments, we make the simplifying assumptions 
that a — 1/2 and v = 0. These assumptions are tantamout to assuming that, under 
uniaxial compression, the material exhibits neither elastic nor plastic expansion in the 
lateral direction. In other words, in this approximation the snow does not interact with 
the container side walls and the confined snow behaves exactly as a free standing column 
where we are dealing with a purely uniaxial stress and strain state. The accuracy of these 
assumptions will be demonstrated in the section on model parameters below. 

With a = l/2 the equivalent stress is given by 


2 ^.dev^.dev . 

CT eq = Vij + “3“ 


(1) 


The plastic strain rate as derived from the associated flow rule is 

e pl = e— (cr dev + ^Tr(<r)lV (2) 

Identifying the compression direction with the x axis direction of a Cartesian coordinate 
system and using v 0, the stress, elastic and plastic strain tensors then have the simple 
structure 

/ 1 0 0 \ / 1 0 0 \ 

er = s I 0 0 0 J = E(e(x))e el (x ), e pl = e(x) I 0 0 0 J . (3) 

\ 0 0 0 / \ 0 0 0 / 

Owing to the parameter choices a — 1/2 and v = 0 there is no lateral expansion or 

contraction of the sample and hence no interaction with the confining walls. The measured 
data for the axial compaction forces and lateral forces allow us to assess the quality of 
this approximation: At the end of each deformation experiment we find typical axial 
forces of 250 N which, given the sample geometry, translates into a typical axial stress 
a xx = 5 x 10 4 Pa. With an end height of the sample of L' « 140 mm, the typical lateral 
force of 45 N translates into a typical lateral stress of a zz = 12 x 10 2 Pa = 0.025a xx . Thus, 
lateral stresses can be neglected to a very good approximation: the elastic state of the 
sample is almost fully characterized by the scalar axial stress s. 
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According to Eq. 3, the volumetric and shear strains are given by £y = —eg 1 = e, hence 
the plastic deformation state is completely characterized by the single scalar variable 
e < 0. The stress equilibrium condition diver = 0 implies that a — const. (x). From this 
condition we can derive a relation between the stress a, the plastic strain e(x), and the 
strain e ext that is externally imposed by the downward displacement of the top surface 
of the sample. Let x be the material coordinate (the position in the sample prior to 
deformation) and let (•) = L~ l J • dx denote the averaging operator on L. We then find 
that 

s — E e g (e ext — (e)), with E~$ — (E~ 1 (e(x))) (4) 

where E e g is the effective elastic modulus of the sample. 

The equivalent stress is cr eq = — s, and the yield condition reads 

f{e, s) = -s- (E R + Eg) exp[—ae] - E ol 2 g^ = 0. (5) 

Eg evolves according to 

d t E s = -Es + —(Eg - E s ). (6) 

£S T 

The system of integrodifferential equations for e, s and Eg is solved as follows: In each 
time step we increase the externally imposed strain according to the corresponding strain 
rate, and evaluate s and e self-consistently from Eqs (4) and (5) with the boundary 
condition d x e(0) = d x e(L) — 0. By comparing with the previous time step we determine 
the local plastic strain rates e, evaluate the changes in Eg from Eq. (6), and repeat this 
process until a given end strain is reached. 


3 Model parameters 

3.1 Elastic constants 

Material parameters of snow are often difficult to measure. This is particularly true 
for elastic properties, where quasi-static measurements are skewed by creep deformation, 
while dynamic measurements may be skewed by pore fluid effects. Furthermore, the 
intrinsic variability of the material leads to a huge scatter of values even if one considers 
a series of measurements with similar methodology and sample preparation [5]. Therefore, 
any values of elastic constants should be considered as indicative. 

Regarding the density dependence of Young’s modulus, we refer to dynamic measure¬ 
ments of Sigrist [6, 7] which cover the density range from p — 200kg/nr 3 to 400kg/m 3 . 
The data can be well represented by a power-law dependenccy, E oc p b with b — 2.94, 
which compares well with the behavior of other solid foams [4], We thus adopt the value 
b — 3 in our calculations. However, the absolute values obtained by Sigrist cannot be 
directly applied to the present context owing to the highly dynamic measurement tech¬ 
nique. We instead refer to measurements of Schweizer and Camponovo [8] who determine 
shear moduli G in a temperature and elastic strain rate range comparable to the present 
experiments. For a series of tests on natural snow with density p 240 kg/m 3 , they give 
a mean shear modulus of G — 0.75 MPa. With Poisson’s number v ss 0 (see below) this 
translates into E ss 2 G = 1.5 MPa and, using the power law to extrapolate, we find for 
snow of density p 370 kg/nr 3 as used in the present experiments a characteristic elastic 
modulus of E ss 5 MPa which we adopt as our initial value of Young’s modulus. 
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There are few measurements of Poisson’s numbers for snow reported in the literature, 
with reported values being typically very small [9]. In the present context we can obtain an 
upper estimate of Poisson’s number from our own lateral force measurements. Identifying 
the side plate normal direction with the z direcation and assuming that the lateral force 
is exclusively due to elastic cross-expansion, we find from cj zz /g xx — 0.025 — u/{l — v) 
that v < 0.025. This allows us to approximate v ~ 0 which yields the simple elastic 
constitutive relation a 7J = Eef,. 

3.2 Parameters of the plasticity model 

Fundamental to our constitutive formulation is the parameter a which governs the vol¬ 
umetric response in the plastic regime. To determine this parameter, we observe that 
the experimental deformation behavior is characterized by the emergence of localized 
deformation in the form of compaction bands that are oriented perpendicular to the 
compression direction. As discussed by [10], such behavior occurs for a — 1/2. We note 
that ol — 1/2 is also quite commonly used in deformation of metallic foams, see the pa- 
rameterizations of [10], [11] and [12]). Again, we can use our lateral force data to derive 
an upper estimate for the deviation 5 — 1/2 — a from this value. For S ^ 0, the lateral 
plastic expansion is, to lowest order in 5 , given by sf z — {2/?,)8e^ x . Assuming that the 
lateral force at the end of the deformation experiment is exclusively due to the elastic 
accomodation of this lateral plastic expansion, and using the elastic modulus value given 
above to evaluate sf z = a zz /E, we find an upper estimate 5 < 10 -3 which shows a 1/2 
to be a very good approximation. 

The parameters characterizing the strain and time dependence of the yield stress in our 
model - the exponent a relating strength to density, the initial strength E = Eg + E R , the 
residual strength S R , the softening strain eg, and the relaxation time r - were obtained 
by matching simulated and experimental stress-strain curves. This yielded the values 
a = 7, E = 5 x 10 4 Pa, E R = 0.25 x 10 4 Pa, es = 0.05, and r = 10 s. The fitted strength 
values are well within the range reported in the literature [5, 13, 14], and the same is true 
for the relaxation time r [15, 16]. 


4 Stability analysis 


For a — 1/2, instability of homogeneous uniaxial deformation with respect to formation of 
localized bands is associated with perturbations of strain in the axial direction, i.e., bands 
are oriented occur at 90 degrees to the stress axis [10]. Instability occurs if the derivative 
of the yield function with respect to the plastic strain becomes positive. Neglecting for 
the time being the gradient term, the stability requirement is thus: 


dj_ 

de 


ds 1 

= —— a(E s + S R ) exp[-ae] -|-E s exp(-ae) < 0. 

Oe £g 


(7) 


Strongly localized perturbations are of negligible influence on the stress given by Eq. 4, 
hence we obtain the simple stability criterion 

£ s 


as s > 


Eg + Ef 


( 8 ) 


This is tantamount to requiring that the rate of hardening due to densification, with 
structural strength fixed, exceeds the rate of softening due to bond breakage, calculated 
with density fixed: Densifcation-induced hardening must outweigh shear softening. 
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Let us first consider the fully aged initial state. In this state, Xs = Xg and (in)stability 
does not depend on the rate of deformation but only on the relative efficiency of hardening 
and softening processes in the fully aged state. For the parameters given above, we can 
see that the initial state of the material is within the unstable domain, hence initial 
deformation is in our simulations always characterized by strain localization. 

We now consider a homogeneous deformation process without spatial dependency 
and fixed imposed strain rate. In this case, a stationary solution of Eq. 6 is Xg ~ 
Xg/(1 + er/es), i.e., strength is reduced due to the on-going destruction of bonds in 
the snow microstructure. This leads to an upper critical strain rate above which the 
specimen behaves after initial softening like a cohesionless aggregate which compacts 
homogeneously. From the parameters above, the critical strain rate can be estimated 
to be of the order of f s -1 , which is well above the strain rate regime accessible in our 
experiments. 


5 Correction for friction forces 

The lateral stresses measured at the end of our experiments are small, which justified 
the assumption of a uniaxial stress state as consistent with the parameters a — 0.5 and 
v 0. Nevertheless, the lateral forces cause friction forces which oppose the movement of 
the anvil and thus influence the measured force data. To clarify the nature of these forces 
we observe that, during passage of the first compaction band, there is an asymmetry of 
the recorded force signal depending on whether the band moves downwards or upwards. 
During motion of a compaction front from top to bottom of the sample, the total force 
rises as more snow becomes compacted and moves with the anvil. By contrast, the force 
remains steady when snow is being compacted but not moving with the anvil (motion of 
a compaction front from bottom to top). Comparing both cases allows us, together with 
the lateral force data, to estimate the coefficient of friction at the wall. 

With 5 — 1/2 — a the lateral plastic expansion of the snow is given by eh' = ef z — 
(2/3 )5e where e is the axial compressive strain evaluated for a — 1/2. Because of the 
confinement due to the side walls, this expansion is offset by a compressive stress a yy — 
a zz — (2/3 )Se(x)E(x) which causes a total lateral force 

y wL(e(x)E(x)) (9) 

If frictional motion occurs at the side walls, i.e., if the displacement rate u x is negative, 
this in turn causes a total friction force 

^^wL(Q(-u x )e(x)E(x)) (10) 

where // is the coefficient of friction at the side wall-snow interface and 0 is Heaviside’s 
function. This force superimposes on the force required for axial compression of the snow 
sample and needs to be taken into account when comparing experimental and simulated 
force distance curves. In particular, during passage of the first band our constitutive 
model predicts a constant stress, hence, any slope of the measured force distance curve 
can be attributed to friction. From the measured slope and the lateral forces, we can 
deduce an estimate for the coefficient of friction at the side walls, which is found to be in 
the range 0.4 < // < 0.6. This is in rough agreement with a direct friction measurement 
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using a weighted block of snow sliding on the side surface of the container, which yielded 
fi « 0.33. 

Our simulated force distance curves shown in the main paper have been corrected for 
the effect of friction. With // ss 0.5, the total friction force at the end of the experiment is 
of the order of 45 N which, while significantly below the total compressive force of about 
250 N, is not negligible. 
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Figure 1: Experimental apparatus (schematically). Snow is enclosed in the lower part 
(a) and compacted by an anvil (b) driven by motor (c). Force is measured by a load cell 
(d); lateral forces can be determined by measuring the bending of the front/back plate 


m (a). 











